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Summary 

The design of a very fast, automatic black-box code for 
homogeneous, gas-phase chemical kinetics problems requires 
an understanding of the physical and numerical sources of 
computational inefficiency. Some major sources reviewed in 
this report are stiffness of the governing ordinary differential 
equations (ODE’s) and its detection, choice of appropriate 
method (i.e,, integration algorithm plus stepsize-control 
strategy), nonphysical initial conditions, and too frequent 
evaluation of thermochemical and kinetic properties. Specific 
techniques are recommended (and some advised against) for 
improving or overcoming the identified problem areas. It is 
argued that, because reactive species increase exponentially 
with time during induction, and all species exhibit asymptotic, 
exponential decay with time during equilibration, exponential- 
fitted integration algorithms are inherently more accurate for 
kinetics modeling than classical, polynomial-interpolant 
methods for the same computational work. But current codes 
using the exponential-fitted method lack the sophisticated 
stepsize-control logic of existing black-box ODE solver codes, 
such as EPISODE and LSODE. The ultimate chemical kinetics 
code does not exist yet, but the general characteristics of such 
a code are becoming apparent. 

Introduction 

The present work is motivated by the need for reliable and 
computationally efficient methods for the numerical modeling 
of continuous combustion phenomena in multidimensional 
reactive flows. 

Prior to about 1970 the literature pertaining to the solution 
of homogeneous, gas-phase combustion kinetics batch reaction 
equations was concerned with single-point calculations 
associated with the determination of chemical kinetic rate data 
by shock tube or flow reactor or with the nozzle performance 
of chemical rocket motors (refs. 1 to 4). In both computational 
scenarios accuracy was the most important concern in selecting 
an integration algorithm. Computational efficiency was of 
concern as well, but evolutionary increases in the execution 
speed of new computers tended to diminish its importance. 
As a result, a number of highly accurate, moderately efficient 
computer codes were developed (refs. 5 and 6), and it could 
be safely said that the single-point kinetics calculation problem 
was solved satisfactorily. 


As attention was focused on modeling of reactive flows in 
the late 1960’s and early 1970’s, it was found that the single- 
point codes were not sufficiently fast, and were unnecessarily 
accurate, for practical application to flows with complex 
reaction mechanisms on large computational grids (refs. 7 and 
8). In particular, whether these flows are being modeled from 
a space-discretized Eulerian approach with finite-rate reaction 
kinetics treated by operator splitting or from a mass-discretized 
Lagrangian approach using the method of fractional steps to 
treat the finite-rate chemistry, there is a common need for a 
moderately accurate, extremely fast, homogeneous batch 
chemistry integrator to give approximate solutions for the 
thousands of resulting initial value problems. The same 
requirement arises in the single-point, stochastic simulation 
of turbulent, inhomogeneous gas-phase continuous-combustion 
systems (refs. 9 to 11). 

Seeking computationally efficient methods for approximately 
solving the stiff, strongly coupled, nonlinear ordinary 
differential equations (ODE’s) governing this problem requires 
that the actual source of the difficulty be recognized, that is, 
whether specific computational problems arise from the 
physics of the system or from an inappropriate choice of 
numerical methods. 

Governing Differential and Algebraic 
Equations 

The system of ordinary differential equations describing 
adiabatic, homogeneous gas-phase chemical reaction at 
constant pressure is given by 


-7* =/,<**, r ) i,k = l,NS (1) 

dt 

where 


fi= ~ P ' 2} ( a ij ~ a u) ( R j ~ R -j> ( 2 ) 

7=1 

where the molar forward and reverse reaction rates per unit 
volume Rj and R__j 9 respectively, are given by 
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R_j = A_jT B -j exp ^ — A / Jl (4) 

In equations (1) to (4) o [ is the mole number of the ith 
species (kmol //kg mixture), NS is the total number of distinct 
chemical species, T is the temperature, and p is the mixture 
mass density. The og’s are the stoichiometric coefficients of 
reactant and product species / in the yth reaction (a-j and 
respectively); J is the total number of independent chemical 
reactions. The quantities A jt A_ Jt B Jt B_ Jt Tj , and T_j are 
constants in the modified Arrhenius expressions, equations (3) 
and (4), for the forward and reverse reaction rates, 
respectively. 

The mixture mass density p in equations (2) to (4) is 
determined by the equation of state for an ideal gas 

P 


where P is the absolute pressure, R is the universal gas 
constant, and o m is the reciprocal mean molar mass of the gas 
mixture, given by 


The purpose of this decomposition is to enable factorization 
of the ith mole number cr, from the destruction term A 

A = Lpi (9) 

where L n the loss coefficient for the ith species, is obtained 
simply by dividing D, by erg 

J 

Li = (pc,)-' £ (otjjRj + a.jR-j) ( 10 ) 

7 = 1 

In this report attention is restricted to constant pressure, 
adiabatic chemical reactions. For such problems conservation 
of thermal energy is expressed by an algebraic equation of 
constant enthalpy as a constraint on equations (1) to (4): 

NS 

^ h,a l = h (} — constant (11) 

;= ! 

where h l is the sensible-plus-chemical molar specific enthalpy 
of the ith species, and h ( , is the mass specific enthalpy of the 
mixture. 

The algebraic equation (11) can be differentiated with respect 
to time so that the enthalpy conservation equation is expressed 
as an additional ODE for the temperature: 



(5b) 


The net production rate of ith species, f in equations (2) 
to (4), may be expressed as a difference between two positive- 
definite terms 


/ = a - a (6) 

where 

J 

Qi = p ] y; ( a ij R ~j + a tj R j) < 7 > 

7=1 

and 

J 

A = P 1 (otijRj T a (S) 

j= l 


(IT 

dt 


NS 

t h ' f ' 

i= 1 

NS 

i= i 


(12) 


where c pi is the constant-pressure molar specific heat capacity 
of the ith species. 

Either equation (11) or (12) can be used in the equation set. 
When equation (11) is used, the temperature is calculated from 
the species mole numbers and the initial mixture enthalpy. An 
iterative technique is employed, and the temperature is adjusted 
until equation (1 1) is satisified. In this method the number of 
independent ODE’s is equal to the number NS of distinct 
chemical species. The use of equation (12) to solve for the 
temperature increases the number of independent ODF/s to 
NS -F 1 . In this method the integrator tracks the solutions for 
both the species mole numbers and the temperature. 

Statement of the Physical Problem 


In equations (6) to (8), Q t and D t represent the gross rates 
of production and destruction of the ith species, respectively, 
because of the contributions of all the J forward and reverse 
reactions in the prescribed mechanism. 


During the process of modeling a chemically complex 
reactive How field, it is typically necessary to solve thousands 
of initial value problems of the following form: given a 
prescribed pressure P, an initial temperature T 0 , and a 
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corresponding initial set of mole numbers ((a/] 0 , i — l, NS), 
find the resulting temperature and mole numbers at the end 
of a prescribed time interval At. 

Figures 1 and 2 illustrate solutions to equations (1) to (12) 
for two typical problems (refs. 12 and 13) of stoichiometric 
combustion of reactive mixtures following rapid compression 
in a shock tube. Three distinct physical and chemical regimes, 
commonly denoted as induction, heat release and equilibration, 
are apparent. 

The induction regime is the period of time immediately 
following some form of homogeneous bulk ignition. In the 
present examples ignition is due to the passage of a shock wave 
through the stoichiometric mixture of fuel and air. During the 
induction period concentrations of reactant precursors and 
intermediate chain carriers (such as O, OH, H, and CH X ) 
increase by many orders of magnitude, from very small initial 
concentrations to values sufficient to initiate endothermic 
reactions. These, in turn, lead to oxidative and thermal 
pyrolysis of the hydrocarbon fuel, which produces CO and 
H 2 . In this regime D { (eq. (6)) is small; f is large and 
dominated by Q r During the early part of the induction 
period, the coupling with the energy equation is weak, so that 
an essentially isothermal reaction is obtained. At the end of 
the induction period, observable ignition occurs, as exhibited 
by an exponentially increasing temperature and accompanying 
rapid depletion of reactant concentrations. 

The induction period ends, and the heat-release period 
begins, when a sharply defined change in temperature and 
molar concentrations occurs. In this regime the full chemical 
mechanism is active, with very strong temperature coupling 
through the enthalpy conservation equation. The heat-release 


period ends after concentrations of the chain-carrying 
intermediates (such as O, OH, H) have reached their peak 
values, and all species have begun to approach their chemically 
equilibrated concentrations. 

The equilibration regime is characterized by the monotonic, 
asymptotic approach of all species concentrations and of the 
temperature toward their chemical equilibrium values. During 
equilibration both Q x and D, are large numbers, but with a 
small difference. The equilibration process does not have a 
clearly defined termination, because of the asymptotic nature 
of the approach to the equilibrium state. However, since the 
equilibrium state can be computed a priori by an efficient Gibbs 
function-minimization scheme (refs. 7 and 14), the end of the 
equilibration period can be defined as the time at which the 
values of all thermochemical variables are within (say) 1 
percent of their equilibrium values. 

Statement of the Computational Problem 

Equations (1) to (12) are to be solved approximately in a 
stepwise fashion by constructing a sequence of approximate 
solutions on the point set h n € [0,A t] starting with the initial 
conditions [a,) 0 , T 0 , and P. The computational mesh [h n ] is 
not given; its construction is part of the task (ref. 15). 

The problem is to find the optimal computational method, 
that is, a numerical integration algorithm and an appropriate 
scheme for controlling local truncation error (accuracy) and 
determining an optimal set of time step intervals \h n }. The 
most efficient method will optimize the tradeoff between 
numbers of steps per interval [0,Ar] and the number of 
computer operations (computational work) per step. 
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Figure 1. — Variation of species mole fraction with time for constant-pressure (10.0 atm) combustion of stoichiometric methane (pyrolized to CO and 
H : ) and air. Initial temperature, 1000 K; reaction mechanism includes 12 reactions and 11 species (ref. 12). 
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Reaction time, s 

Figure 2.— Variation of species mole fraction with time for constant-pressure 
(2.0 atm) combustion of stoichiometric hydrogen and air. Initial temperature, 
1500 K; reaction mechanism includes 30 reactions and 15 species (ref. 13). 


In principle, classical explicit integration algorithms, such 
as the Runge-Kutta family, may work satisfactorily on some 
part of the computational mesh. However, implicit algorithms 
require less computational work for the same accuracy and 
so are preferred. Explicit algorithms are not used at all, except 
as predictors in a predictor-corrector formulation. 

Implicit algorithms for stepwise approximate integration of 
equation (1) may be written in component form (with the 
species index subscript suppressed for clarity) as follows 
(ref. 15): 

",H 1 = h,< . ift, + !./„) I + in 

where n denotes the time level at which the approximate 
solution is known, n + 1 denotes the advance time level at 
which the solution a„ H is sought, h n + { is the step length 
(/„+ j - f„), f n+[ is the net production rate at t n + Jf f i and 
\p n are constants or algebraic coefficients characteristic ot the 


particular algorithm chosen, and \J/ n contains all the 
previously computed information from time level t n . In the 
notation of equation (13) an explicit algorithm would be 
represented trivially as a n + 1 = i p n . 

Some of the choices that have to be made in selecting the 
best implicit methods are whether to converge the corrector 
equations by successive substitution (Jacobi, Gauss-Seidel, or 
Jacobi-Newton iteration (ref. 16)) or by simultaneous 
substitution (Newton-Raphson iteration), and if by Newton- 
Raphson iteration, whether to choose pivotal Gauss 
elimination, LU-deeomposition (ref. 17), or Hessenherg 
decomposition (refs. 18 and 19) to carry out the work-intensive 
matrix operations required. The local truncation error must 
he estimated in order to assure acceptable accuracy, which 
can be done by the full-step, half-step method, from the 
difference between predictor and corrector solutions, or from 
the difference between two different-order correctors 
(ref. 20). If step size is to be controlled by both accuracy and 
iterative convergence efficiency, the convergence rate must 
also be monitored. All these choices are interdependent and 
highly problem-dependent, so that a single, best, all-purpose 
computational method for all problems of the form ot 
equation (1) does not exist yet. 


Stiffness Problem 

An efficient, automatic computational method must be able 
to identify stiffness of the governing system of differential 
equations and then be able to invoke the appropriate numerical 
integration algorithm (refs. 18, 19 and 21). Unfortunately 
numerical and computational specialists do not agree on a 
single technical definition of the term “stiff,” as applied to 
systems of linear or nonlinear ODE’s (ref. 22). However, 
Lambert’s (ref. 15) and Shampine’s (ref. 22) definitions are 
pragmatic, self-consistent, and understandable: 

( 1 ) A system of ODE’s is stiff if the exact solutions are stable 
(for a specific set of initial conditions) in the forward direction, 
but the exact solutions in the reverse direction (decreasing time 
or axial direction) are unstable. Any physical system which 
has a strong directional drive toward some form of equilibrium 
state is described by stiff ODE’s. Typically the exact solutions 
of a stiff system of ODE’s exhibit the general appearance of 
decaying exponential functions. 

(2) A problem (defined as the ODE’s, a set of initial 
conditions, and a stepwise numerical method for their 
approximate solution) is computationally stiff if in the region 
of interest the computational step size must be reduced severely 
from that value which would yield the required accuracy. 

fhe most important and useful distinction between the two 
definitions is that stiffness is an inherent attribute of the 
physical system of interest (and, by definition, of its associated 
system of ODE's), whereas computational stiffness is just a 
symptom of having chosen an inappropriate or suboptimal 
scheme for the approximate numerical solution of the ODE’s. 
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If the physical system of interest does not exhibit strongly 
directional, stable behavior— that is, if its exact solutions are 
weakly stable or unstable— then the system and its describing 
ODE’s are nonstiff. Systems with oscillatory exact solutions 
may alternate between regimes of stiff and nonstiff behavior. 

It can be seen by inspection of figures 1 and 2 that during 
equilibration the ODE’s are stable, and the chemical species 
and temperature asymptotically approach their equilibrium 
states. Therefore, using the definitions of stiffness, we can 
classify the equilibration regime as stiff. During induction, 
however, many species and the temperature have positive time 
constants, which indicate unstable ODE’s. The induction 
regime is, therefore, nonstiff. Since the heat-release regime 
is not clearly stiff, it must be regarded as nonstiff as well. 

Nonstiff Integration Algorithms 

In the nonstiff regime, the differential equations are unstable 
(or at best weakly stable), so that small step sizes are required 
regardless of the accuracy of the integration algorithm. 
Therefore, since a limited convergence radius is not the 
restricting consideration, the computationally inexpensive, 
point-iteration methods are preferred to matrix iterative 
methods. 

With the decomposition of f into gross production and 
destruction terms of equation (6), Jacobi-Newton iteration 
(ref. 16) should be used to achieve the fastest convergence 
for equation (13): 

^0.(5+^) __ ^n+lfffl+l/n+l — Pff+l + $ n 

1 3 “ hn+ \&n + \L>n + 1 

and 

= a‘!>, + Aa< s + V> (15) 

where Act is the iterative correction, and ( 5 ) denotes the 
iteration counter. Equations (14) and (15) are iterated until 
converged, that is, until the absolute value of Act is less than 
some suitably small criterion e. 

The loss coefficients L, in equation (14) approximate the 
diagonal terms of the full Jacobian matrix, equation (17). They 
are used in equation (14) only to accelerate convergence, so 
that they need not be reevaluated at each iteration. 

The step size h n+l is controlled so that the rate of 
convergence (defined as Act^/VAct^j) is about 0.5, which, 
experience shows, optimizes the number of steps in the 
prescribed interval and the number of iterations per step 
(refs. 18 and 19). Gear's nonstiff method (ref. 23) can be 
employed to select automatically the lowest order Adams 
predictor-corrector formula which will satisfy a prescribed 
level of accuracy, as determined by an estimate of the local 
truncation error. The packaged ODE solver code EPISODE 
(ref. 24), with method flag MF = 13, incorporates the nonstiff 
method recommended. 


Automatic Detection of Stiffness 

If the nonstiff integration algorithm cannot achieve the 
desired rate of iterative convergence, the problem has become 
computationally stiff, and a change to a stiff integration 
algorithm must be considered in order to maintain optimal 
computational efficiency (refs. 21 and 22). However, it is not 
clear in practice where the optimal changeover point occurs. 
In the context of the present problem, a reliable criterion is 
available from the problem physics: when two or more species 
are in quasi-steady-state (defined for this purpose as occurring 
when the absolute difference between Q { and D { of eqs. (6) 
to (8) is less than one thousand times their sum), the system 
is considered to be computationally stiff, and the use of a stiff 
integration algorithm is indicated. 

Stiff Integration Algorithms 

In stiff regimes, the large step sizes admitted by accuracy 
requirements are too great for convergence of the best point- 
iterative schemes. Therefore, it is necessary to use some form 
of Newton iteration to converge equation (13), which in turn 
requires the evaluation of an exact or approximate Jacobian 
matrix and the iterative solution of matrix systems of 
dimensions NS by NS\ the latter requires an enormous increase 
in computational work per step (ref. 15). 

A suitably modified form of equation (13) for computing 
the Newton-iterative corrections may be written in vector form 
as 

(/-/w./wfya iVV* 

= A» + l& + |/ftl -?Wl +*„ ( 16 ) 

where / is the diagonal unit matrix and J is the Jacobian matrix, 
defined in component form by 

j,k = = (p^r 1 E (“0 ~ a ‘j) ( R j a ki - R -j a kj) (1 7 ) 

° k j = 1 

and where the corrections Act, are to be used as before (see 
the Section Nonstiff Integration Algorithms). In equation (17) 
the partial derivatives with respect to the density are assumed 
negligible in comparison with the other terms. 

A key point in the solution strategy is that the notation J 
for the Jacobian matrix in equation (16) purposely omits 
reference to either the time level or the iteration counter. The 
idea is to use an old Jacobian as long as possible, to avoid 
recalculation of J by equation (17), and more significantly, 
to avoid repeated decomposition of the /-matrix. The iterative 
convergence rate is monitored as in the nonstiff case, but 
inefficient convergence rates (>0.5) are taken as a signal to 
reevaluate J and to decompose the new matrix. When this 
iterative convergence is used in connection with the family 
of variable-order, backward difference integration algorithms 
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(Gear’s stiff method (ref. 23)), a very efficient stiff method 
results. 


Evaluation of Thermochemical and 
Kinetic Rate Data 


Nonphysical Initial Conditions 

In a single-point calculation, such as the shock-initiated 
combustion problems illustrated in figures 1 and 2, the initial 
mole numbers of reaction intermediates and of products are 
usually assumed to be equal to zero, rather than the very small 
but nonzero initial values which actually exist. In 
multidimensional reacting flow calculations, initial mole 
numbers are determined by a weighted averaging of mole 
numbers over adjacent grid nodes or from spatially adjacent 
fluid elements, depending on whether an Eulerian or a 
Lagrangian description of the flow is used. Since the properties 
associated with some neighboring nodes or fluid elements may 
represent chemically and thermally hot (ignited and burning) 
states, while others may represent cold (unignited) reactant 
states, the warm mixture properties resulting from the weighted 
averaging may be physically meaningless. 

When this happens, a numerical imbalance occurs, for one 
or more species, between the production and destruction terms 
in equation (6), (3, and D n respectively. A spuriously high 
positive or negative species net production rate /, : is the result. 
The effect on any numerical integration algorithm is that 
unnaturally small step sizes are required to resolve the very 
large predicted change in mole numbers due to the 
nonphysically high species reaction rates, which results in 
excessive computational work. 

The solution is to “filter” the initial conditions in order to 
provide physically meaningful initial mole numbers and net 
species production rates. One method that accomplishes this 
filtering is as follows: 

(1) Choose an initial (i.e., first) time step /Z| = [max(L,)| '. 

(2) Calculate a predictor set of mole numbers by using the 
explicit “filtered Euler” approximation (ref. 25), 


(7 


(0) 

/.I 




1 - exp( - L i () h i) 
L U) h i 


(18) 


where the superscript (0) indicates the result of the predictor 
step 

(3) Iterate an implicit Euler corrector to convergence: 

' ) — < y /.o + h\ 7/.V (19) 


A preprocessor such as CREK (ref. 26) or CHEMKIN 
(ref. 27) is used to read the kinetic mechanism, the kinetic 
rate data constants required in equation (3), and the fifth-order 
polynomial fit coefficients required to evaluate the 
thermochemical properties required for the calculations in 
equations (2) to (12). If the reverse rate data A_ jt B r and 
T_j of equation (4) are not prescribed, the reverse rate 
constants must be calculated from the following detailed- 
balance relation: 


ss 

, , v (“*/ ■■<%) 

k_i - k ; (Riy 1 exp 


NS 

E («.v 


/= I 


a ^RT 


(20) 


This calculation requires repeated evaluation of the molar 
specific enthalpy //, and the molar specific 1-atm entropy s]' 
in order to evaluate the /th species 1-atm specific Gibbs 
function g," 

H? = hi- Ts? (21) 


If the loss in accuracy is acceptable, the reverse rate 
constants can be precalculated over a range of expected 
temperatures and obtained by a least-squares fit to the reverse 
rate constants as given by equation (20) (ref. 26). 

Another important technique is to avoid unnecessary 
evaluation of thermochemical and rate data by locally 
linearizing the rate constants and thermochemical data, so that, 
during the course of iterative convergence of equation ( 14) 
or (17), the thermochemical and kinetic rate data are not 
reevaluated while the current temperature is within a local 
window (T, T + | AT| ). The size of this window is controlled 
so that resulting errors in the approximate solution observe 
the prescribed error bound e on local truncation error in the 
following way: 


AT ~ - 


max 

T. 

+ B, 

./ 

7 1 


cT 

T_ j 


+ B 


( 22 ) 


Use of this control strategy has been shown to reduce the 
total computational work significantly (refs. 13 and 28). 


The predictor equation (18) damps the bad initial rates, and 
the corrector equation (19), being fully implicit, has no 
memory of the initial rates at all. The filtered set of mole 
numbers o l is now ready for the main integration routine. 


Prospects for Further Improvements 

No code presently exists which combines all of the speed- 
enhancing techniques discussed in preceding sections. Because 
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of the often subtle interactions between parts of an integrated 
code, a certain amount of trial-and-error combination of 
techniques seems to be an inevitable requirement for further 
improvement. In this section prospects for further speed 
reductions in each of the three physical regimes— induction, 
heat-release, and equilibration— are considered separately, and 
finally a development path toward the ultimate integrated code 
is suggested. 

Induction Regime 

The increasing-exponential behavior of reactive species 
apparent in figures 1 and 2 suggests the use of exponential 
functions (or their rational-function or Pade approximants 
(ref. 29)) as interpolants, rather than the classical polynomial- 
interpolant methods, in determining the /3’s and \t’s in equation 
(13). The exponential-fitted algorithms (refs. 29 to 32) do, 
in fact, interpolate the reactive species with more accuracy 
for the same computational work than the Adams methods; 
however, the lack of an efficient step-size control strategy 
cancels much of the advantage (refs. 13, 28, 31, and 32). 

Other methods have been suggested to save computational 
work in the induction regime, including formally replacing 
time with temperature as the independent variable and directly 
calculating the induction time (ref. 12) or formally replacing 
time with the mole number a l of one of the reactive species 
as the independent variable and integrating the governing 
equations in the phase plane. Neither of these two proposals 
has been demonstrated to be computationally efficient. 

Heat-Release Regime 

There does not appear to be any obvious way to increase 
significantly the speed of computation in this intermediate 
regime, which is neither clearly stiff nor nonstiff. The 
exponential-fitted methods offer no apparent advantage over 
polynomial-interpolant algorithms in this regime, because of 
the pathological variations of species mole numbers with time 
(extrema and inflections). Perhaps a high-order, single-step 
Obreschkoff (spline) method (ref. 33) would give performance 
superior to that of the typically multistep, low-order stiff or 
nonstiff algorithms featured in Gear’s methods (refs. 23, 24, 
and 34). The implicit or semi-implicit Runge-Kutta methods, 
known also as Rosenbrock methods, may offer some special 
advantage in this regime (refs. 35 and 36). 

Equilibration Regime 

Because exponential functions with negative time constants 
exhibit asymptotic-decay behavior, and even high-order 
polynomials do not, the exponential-interpolant algorithms 
appear to be inherently more accurate than polynomial- 
interpolant algorithms of comparable computational work per 
step. Brandon (refs. 31 and 32) claims approximately sixth- 


to eighth-order accuracy for the same work as that of a second- 
order Adams or backward-difference algorithm. The 
developmental batch kinetics code CREK1D (refs. 37 and 38), 
which uses the exponential-fitted trapezoidal rule integration 
algorithm (refs. 29 to 32), has been demonstrated to give 
comparable performance for accuracy equal to that of the best 
Gear’s method code available, LSODE (refs. 13, 28, and 39). 
However, because of the inefficient step-size control strategy 
currently used in CREK1D, LSODE is, at present, the better 
method of the two for use in the equilibration regime. 

Other methods proposed to reduce the computational work 
associated with the matrix operations required for Newton 
iteration include using perturbed functional iteration (ref. 40) 
and reducing the size of the computational matrix by using 
the kinetic rate data for only three-body reactions as constraints 
on the sum of the mole numbers a m in a constrained 
equilibrium method (ref. 41). The far more restrictive, but 
time-honored, assumption of quasi-steady-state (QSS) behavior 
for reactive species is not recommended, as the time saved 
by a modest reduction in size of the matrix is largely offset 
by the difficulties in solving the mixed, differential-algebraic 
system of equations; in any case, the QSS assumption is not 
element-preserving, so that periodic adjustments are required 
to restore conservation of elemental mass (ref. 42). 

Toward the Ultimate Code 

The ultimate batch kinetics code, when called from a reactive 
flow field modeling code, would automatically (i.e., in a black- 
box manner) perform the following actions: 

(1) Filter the initial conditions 

(2) Identify three physical and computational regimes: (a) 
induction (unstable nonstiff), (b) heat-release (mixed stiff and 
nonstiff), and (c) equilibration (stable stiff) 

(3) Cali the appropriate integration method (algorithm plus 
step-size control logic) 

(4) Reevaluate thermochemical and kinetic rate data only 
as often as necessary to observe prescribed error bounds 

Ideally more than one integration method would be available 
for each of the three regimes, depending on whether the user 
elects a high, medium, or very low accuracy or computational 
work. At the present time LSODE (ref. 34) offers the best 
available method for all three regimes, but there is clearly room 
for further improvement. 

Concluding Remarks 

The design of a very fast, automatic integration code for 
homogeneous, gas-phase chemical kinetics depends on 
understanding the physical and numerical sources of 
computational inefficiency. Some specific techniques for 
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overcoming three of the major sources of inefficiency- 
stiffness and computational stiffness, nonphysical initial 
conditions, and unnecessary evaluation of thermochemical and 
kinetic properties— have been recommended. 

The ultimate black-box code will automatically filter the 
initial conditions, select the best integration method for the 
physical and computational regime identified, and avoid 
unnecessary calculation of thermochemical properties. 

Progress remains to be made in the areas of further reducing 
or eliminating altogether the burdensome work of matrix 
calculations and in devising very low accuracy or approximate 
integration methods. However, with the sources of present 
inefficiencies now reasonably well understood, it is clear that 
further improvements in computational efficiency are possible. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, February 11, 1986 
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